The effect of network structure on phase transitions in queuing networks 
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Recently, De Martino et al 0, @] have presented a general framework for the study of transporta- 
tion phenomena on complex networks. One of their most significant achievements was a deeper 
understanding of the phase transition from the uncongested to the congested phase at a critical 
traffic load. In this paper, we also study phase transition in transportation networks using a dis- 
crete time random walk model. Our aim is to establish a direct connection between the structure 
of the graph and the value of the critical traffic load. Applying spectral graph theory, we show that 
the original results of De Martino et al showing that the critical loading depends only on the degree 
sequence of the graph — suggesting that different graphs with the same degree sequence have the 
same critical loading if all other circumstances are fixed — is valid only if the graph is dense enough. 
For sparse graphs, higher order corrections, related to the local structure of the network, appear. 
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I. INTRODUCTION 

During the past few decades, the physics community 
has witnessed enormous progress in the research on com- 
plex networks 0, H[. Transport processes on networks 
represent an important class of dynamical systems with 
a wide range of applications, including data traffic on 
the Internet, vehicle traffic on highways, virus spread 
between hosts, or rumour spread in social networks, to 
name but a few. A simple, general model can be used 
to describe these dynamical systems quite accurately, in 
which particles are transported between the nodes of a 
network. 

In certain transport networks the particles are served 
by queues, residing at the nodes of the network. One 
of the most prominent examples is the Internet, where 
the data packets play the role of the particles. Queuing 
networks exhibit several interesting phenomena, for ex- 
ample, below a critical traffic intensity the system is in a 
free state and the average queue length fluctuates around 
a finite value. Above a critical traffic load, however, one 
or more queues become congested and the average queue 
length diverges. 

Various models have been developed to model the In- 
ternet traffic. Deterministic and probabilistic routing 
strategies were compared in Q . The authors of @ stud- 
ied a shortest-path routing model where the probability 
of packet transmission depended on the queue lengths. 
The authors of 0] studied how traffic congestion is af- 
fected by the capacity of the nodes in a simple shortest 
path routing model. A more elaborate routing strategy 
was studied in [8|, where the packets were forwarded to 
the neighbor that minimized an effective distance to the 
packet's destination. Beyond that, several attempts have 
been made to find routing strategies that are less sensi- 
tive to congestion f9rJlTl|. 

Packet level simulations of these models clearly indi- 
cate phase transitions between the free and congested 
phases, and several characteristics of the queuing net- 



works exhibit power law dependence from the traffic in- 
tensity close to the transition point 0-0, E2 ■ The an- 
alytic description of these models, however, is rather 
limited, because even the simplest routing mechanism, 
namely the shortest-path routing, introduces non-local 
transport dynamics to the system. 

In recent years, extensive research has been undertaken 
to discover the relationship between the topological prop- 
erties of networks and the behavior of the dynamical pro- 
cesses on them [l3l - |l5j . A fundamental question is if 
a dynamical system shows phase transition phenomena 
how does the structure of the network affect the phase 
transition. In models with non-local transport dynamics, 
however, the analytic description of the phase transition 
has only been established for a few special networks, for 
example lattices [!, [l6[ or Cayley trees @, 0, E2 • 

In recent papers by De Martino et al [3,[2|, the authors 
studied the congestion phenomena in arbitrary networks. 
The authors introduced a traffic aware congestion control 
mechanism in their model, and modelled the transport 
process with a simple random walk process, instead of a 
non-local routing mechanism. It has been shown that the 
model exhibits both first and second order phase transi- 
tion, depending on the parameters of a congestion control 
mechanism. 

The main focus of our paper is to gain a deeper un- 
derstanding of the relationship between the structural 
properties of the underlying graph and the congestion 
phenomena, that is the dynamics, and to give a generic 
description of the phase transition point in an arbitrary 
network. Our model is similar to the one presented by 
De Martino et al [1, 2]. We approximate the particle 
transport process with a discrete time random walk pro- 
cess, where packets are generated, absorbed and move 
randomly in the network. Moreover, we assume that the 
delivery of the particles is locally homogeneous, that is 
the probability that a particle will be delivered from any 
node to its neighbour is uniform, and we also assume that 
the queues are independent 
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Our work differs from [l|, |2fl in two important aspects, 
however. Firstly, we do not consider any traffic control 
mechanism in our study. Secondly, instead of using time 
evolution equations, we apply a mean field approximation 
in the long time limit and use spectral graph theory fl7| 
to connect the structure of the network with the traffic 
dynamics. 

The paper is organized as follows. After presenting our 
model in detail in Sec. UH we derive relationships that 
connect the critical traffic loading with the parameters 
of the model (Sec. IIII[) . In Sec. II VI and Sec.[V] we discuss 
the theoretical findings and we compare them to particle 
level numerical simulations and numerical computations. 
Finally, we conclude our work in Sec. IVII Some of the 
details of the analytic calculation are presented in the 
Appendix. 

II. THE MODEL 

We model the transportation network by a simple con- 
nected graph with N nodes and M edges. Moreover, the 
dynamics of the transport networks are modeled by a dis- 
crete time stochastic process. The rules of the stochastic 
process are the following (see Fig.[I|. At each node of the 
graph, there is a queue with infinite buffer capacity. In 
each time step, the first particle of each non-empty queue 
leaves the queue. A particle at node i will be either ab- 
sorbed (i.e. leaves the queuing network) with probability 
fj,i ) or it is delivered to another queue at node j, adjacent 
to node i, with probability Pji. 

The probability that a particle will be absorbed at node 
i can be expressed by /Xj = 1 — . Pji . We will assume 
that the transition probability is constant in time. In 
locally homogeneous network dynamics, the transition 
probabilities can be given by 

p * = if' (1) 

where di is the degree of node i. 

Note that as long as only the queue length statistics 
are concerned and not the fate of individual particles (e.g. 
trajectories or travel times), the order in which the parti- 
cles leave the queues is irrelevant. Therefore, individual 
particles can be considered to be indistinguishable, and 
not only the first, but any packet can be selected from 
the queues for delivery. 

In each time step, after the delivery or absorption of 
the existing particles in the system, new particles can 
also enter the queues randomly. We assume that the 
probability, pi, that a new particle enters the system at 
node i is also constant in time. In addition, we will as- 
sume that the queuing system is open, that is particles 
are generated and absorbed with non-zero probability. 

The queuing network can be in either a free or a con- 
gested state. The network is in free state if, after a tran- 
sient period, the number of particles in the system fluctu- 
ates around an average value. This stationary behavior 




FIG. 1. A schematic figure illustrating the dynamics of the 
model. Particles are generated and absorbed at node i by 
probability pi and fu, respectively. The probability that a 
particle is delivered from node i to node j is denoted by Pji. 
(Color online.) 



does not depend on the initial distribution of the length 
of the queues (l8l.[l9|. In this case, the average number of 
particles arriving to the system equals the average num- 
ber of particles leaving the network. On the other hand, 
in the congested state the average number of particles ar- 
riving to the system is greater than the average number 
of particles absorbed. Therefore, in the congested state, 
the average number of particles in the network will al- 
most surely increase in time. This observation suggests 
the definition of the order parameter 

, ,. n(t+l)-n(t) 
V(Pi,---,Pn) = hm — , (2) 



which measures the expected growth rate of the number 
of particles in the system, n(t), at time t, relative to the 
arrival rate of incoming particles, Y^,iPi B 

In the case of a stationary state, the order parame- 
ter is obviously zero, whereas in the congested state it 
is greater than zero. The transition between the free 
and congested states can be characterized by the critical 
probability, p c = (pi,p2, ■ ■ ■ ,pn) T , where the expected 
arrival rate of the incoming particles equals the expected 
rate of absorbed particles at least at one of the queues. 
It has been shown earlier [2J, |25| that several character- 
istics of these networks (e.g. probability distribution of 
delay times, queue length distribution, etc.) show power 
law dependence on the loading probability near the criti- 
cal point, which suggests a close analogy with the theory 
of phase transitions. 
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III. CONGESTION IN ARBITRARY 
NETWORKS 

In this section we present an analytical estimation of 
the critical point using a mean field approximation 

Let us suppose first that the queuing system is in equi- 
librium. In this case, the expected number of particles 
arriving at each queue is equal to the expected number 
of particles leaving the queue, that is 

&=Pi + A<, (3) 

where £j denotes the expected number of particles leaving 
the queue, A, denotes the expected number of particles 
arriving to node i from its neighbors in one time step, 
and pi is the arrival rate of the particles at node i. Since 
either zero or one particle leaves the queue in each time 
step, £i is also equal to the probability that the queue 
of node i is not empty. Therefore, in the mean field 
approximation, where the queue lengths are independent, 
Aj can be calculated as 

3 

With standard vector and matrix notations we obtain 
that in the uncongested, stationary phase the state vector 
£ has to satisfy the equation (E — P) £ = p, where E 
denotes the identity matrix. Since the matrix E — P 
is invertible (see Appendix A), the loading probabilities 
uniquely determine the components of the state vector £ 
in the uncongested phase: 

^ = (E-P)- 1 p- (5) 

The state of the network can be classified according 
to the state vector If the components of £ satisfy 
the inequality £j < 1 for all i, then the network is in 
an uncongested state, whereas if there is at least one 
node where £j = 1, the network is in the congested state. 
Therefore, based on this condition, the order parameter 
of the system can be calculated theoretically. 

Note first that the balance equation ((3j) cannot hold at 
the congested nodes of the network, because the expected 
number of incoming particles is greater than one, which 
is the maximum of the expected number of outgoing par- 
ticles at a queue. Therefore, the congested queues grow 
steadily, and these queues are never empty. It follows 
that £i = 1 for the congested queues, and the expected 
growth rate of these queues can be given by pi + Aj — 1. 
Based on these observations we can develop an algorithm, 
presented in Appendix[C] that can be used to calculat the 
order parameter numerically for arbitrary networks and 
traffic load. 

In order to validate our model, we compared the order 
parameter obtained from our algorithm with packet level 
simulations on the same network. For the comparison we 
used Barabasi-Albert (BA) (21, Erdos-Renyi (ER) [27| 
and Watts-Strogatz (WS) (28| networks. The loading 
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(b) ER network (per = 0.4) 
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(c) WS network (z = 16 and q ws = 0.1024) 

FIG. 2. The order parameter as a function of homogeneous 
loading probabilities. The graphs show numerical calcula- 
tions (solid line) and particle based simulations (black points) 
on three distinct networks. All graphs had TV = 500 nodes. 
(Color online.) 



probability was the same at every node of the network, 
and both the absorption probabilities and the elements 
of the transition matrix P were random numbers dis- 
tributed uniformly between zero and one. 

The results are shown in Fig. [5J It can be seen that the 
theoretical curve, computed by our numerical method, 
fits very well to the values of the order parameter deter- 
mined by simulations. 

We also validated our results for inhomogeneous traffic 
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FIG. 3. Comparison of the order parameter computed numer- 
ically with our algorithm and obtained from simulation with 
inhomogeneous loading probabilities on an ER graph(iV = 
500 and per, = 0.4). (Color online.) 

loads. For this purpose we generated several realizations 
of random loading vectors, p, and transition matrix P on 
the same ER graph and compared the order parameter 
calculated by simulations and numerical computations. 
Results are shown in Fig. [3] The simulations agree with 
our numerical method very well. 

IV. CRITICAL TRAFFIC LOAD IN 
ARBITRARY NETWORKS 

The main difficulty of using Eq. ([5]) is the computa- 
tional complexity of inverting the matrix E P. More- 
over, even if the matrix can be inverted numerically for a 
particular network, the dependence of the critical point 
on the network structure and the traffic load remains 
obscure. In the case of large irregular networks approxi- 
mations are needed to describe the phase transition an- 
alytically. 

For a detailed analysis of topological effects on the 
critical load p c , let us consider networks with random- 
walk-like particle transport with homogeneous absorp- 
tion probabilities. In this case, if nodes i and j are con- 
nected in the network, the transition probability from 
node i to node j is Pji = (1 — /u)/d,, where d{ is the de- 
gree of node i, and /i is the absorption probability. Using 
the adjacency and degree matrix of the network, P can 
be written in a very compact form: 

P = (1-M)AD _1 . (6) 

Using the spectral decomposition of the transition ma- 
trix, we can write 

N 

AD 1 =Y^^v k w T k . (7) 

k=l 

where Wk and Vk are the left and right eigenvectors of 
AD -1 , respectively, and is the corresponding eigen- 
value. Note that the symmetric matrix D _1 / 2 AD -1 / 2 



has the same eigenvalue spectrum as the transition ma- 
trix AD -1 . Indeed, if Uk is an eigenvector of the later 

1/2 

matrix with eigenvalue Xk, then Vk = D ' Uk will be the 

1/9 _ 

right, and Wk = D ' Uk will be the left eigenvector of 
AD" 1 with the same eigenvalue 1291. 

The transition matrix AD -1 is well known in the the- 
ory of random walks on graphs 1291. Its eigenvalues sat- 
isfy the inequality — 1 < Afc < 1 [29j. The largest eigen- 
value is always Ai = 1 and it is degenerate only if the 
graph is not connected. Moreover, the smallest eigen- 
value is equal to —1 if and only if the graph is bipartite 
[2{|. In the following we will assume that the graphs 
under study are connected. 

It is easy to see that the normalized eigenvector 
of D _1 ^ 2 AD -1 ^ 2 corresponding to Ai = 1 is U\ = 
D 1 ' 2 l/v / 2~M, where M is the number of links in the net- 
work, and 1 = (1, 1, . . . , 1) T . It follows that the left and 
right eigenvectors of AD 1 corresponding to Ai = 1 are 
w x = 1 T /V%M and v 1 = d/V^M = Dl/\/2M, respec- 
tively. 

Using the power series expansion of the function 1/(1 — 
x) and the spectral decomposition of the transition ma- 
trix, we can formally calculate the inverse of E — P: 

N 1 

(E - P)- 1 = J2 r^i TT^wl. (8) 

^ 1 - (1 - M )A fe 

The above sum can be split into three parts that have 
behave differently as /i is varies: 

E T 1 ^ A 1 rp 
vi.wi H h > ; r — Vkwi. (9) 
k h u2M ^ 1 - (1 - u)X k k y ' 

The first term, corresponding to the eigenvalues A^ = 0, 
is independent of /1. The second term, which belongs 
Ai = 1, becomes singular as fi —¥ 0. The last term, 
which consists of all the eigenvalues that are neither zero 
nor one, is finite for every values of \x. 

V. DISCUSSION 
A. Low absorption levels 

In the case of low absorption levels, the second term 
dominates in Q. Therefore, we obtain 

e = (E-p)- 1 pei ^, (10) 

d M 

where p = ^iPi/N and d = 2M/N. Moreover, if the 
loading is homogeneous, i.e. p — pi, the critical loading 
probability is 

d 

Pc = -j (11) 

"max 

where d max is the maximal degree in the graph. 
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Note that in the low absorption limit the critical point 
depends only on the relative spread of the degree se- 
quence, i.e. d max /d, and not on the absolute scale of 
the degrees. In particular, in the case of regular graphs, 
where each node has the same degree, the relative spread 
is d max /d = 1, so the critical point, p c = /i, is indepen- 
dent of the degree of a regular graph. Furthermore, d max 
is always greater than or equal to d, and equality holds 
iff the graph is regular. Therefore, the critical point p c is 
the highest in regular graphs at a given absorption level. 



B. High absorption levels 

If the absorption probability /i is close to one, then 
(1 — /i)Afc is close to zero, so the denominators of the 
third term in @ can be approximated by one. Since the 
eigensystem of AD 1 is complete, we obtain that 



(E-P) 



E 



1 - n dl 1 



H 2M 

and the state vector can be approximated by 

€^P+(i-# 

d H 



(12) 



(13) 



Note that we kept the l//x factor in (fl"3"|) . This way 
Eq. (p~3|) can reproduce Eq. (|10|) in the small absorption 
limit, since p can be neglected compared to the second 
term, which diverges if /i — > 0. 

In the case of homogeneous loading the critical loading 
probability is 



Pc 



/' 



/U + (1 - yu)rf max /d 



(14) 



which depends again only on d max /d, the relative spread 
of the degree sequence. It is remarkable that the critical 
point is almost completely independent of the fine de- 
tails of the network structure in both the small and high 
absorption limit. 



C. Intermediate absorption levels 

Numerical simulations presented later in this section 
show that Eq. (flU)) is valid not only for small and large 
values of fj, 7 but also for intermediate values if the edge 
density of the graph is large. The reasons are the fol- 
lowing. The approximation that leads to Eq. (flU)) is the 
assumption that the term (1 — /i)Afc is close to zero for 
those Afc, that are neither equal to zero nor one. This 
approximation is valid not only when /i is close to one, 
but also if the second largest absolute eigenvalue is close 
to zero. 

If the graph is non-bipartite, the second largest ab- 
solute eigenvalue is related to r, the characteristic time 
until a particle reaches the stationary distribution of a 



random walk, as 1/r = max{|A2|, |Ajv|}. Numerical cal- 
culations, presented in Section IYD1 show that 1/r de- 
creases when the average degree increases and remains 
constant with small fluctuations if the average degree is 
fixed. Therefore, in the case of dense networks, it is plau- 
sible to use the lowest terms of the power series expansion 
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l-(l-/x)A fe 



i + (i-p)\ k + (i-fiyxi + 



(15) 



in the third term of Eq. ([9j). The zeroth order term of 
the series expansion reproduces Eq. (|13|) . If the network 
is not dense, like many real networks, or the absorption 
level is in the intermediate range we need to consider the 
first order term in the power series expansion (I15[) as well. 
Since bipartite and non-bipartite graphs are qualitatively 
different, we will discuss the two cases separately. 



1. Non-bipartite graphs 

It is easy to see that the first order correction to 
(E - P) _1 is (1 - /I )(AD" 1 - dl T /2M). Using this cor- 
rection we obtain that the state vector can approximated 
by 



^p+(l^)AD-V+(l-ri 2 ^. 

d M 



(16) 



In order to study the effects of the topology on the 
critical traffic load, let us consider homogeneous load- 
ing probabilities. After straightforward calculations we 
obtain that in this case the ith component of the state 
vector is 



V 



(1 



- fi) 2 d % di 
M d hi 



(17) 



where hi denotes the harmonic mean of the degree of the 
neighbors of the ith node, 



(18) 



and Af(i) is the set of neighbors of node i. Consequently, 
we obtain that the critical loading probability is 




1 

Pc 



1 



(1 - M) 2 di 
/U d 



+ (l-/i) 



(19) 



which is one of the main results of our paper. 

The main novelty of Eq. (fl"9)) is that it shows that the 
critical point of the phase transition is determined not 
only by the spread of the degree sequence, d max /d, like 
in the low and high absorption limit, but also by hi, which 
depends on the local structure of the network. 

We validated our result on BA and ER networks. We 
calculated the critical load p c by inverting Eq. ([S]) nu- 
merically and then compared the result with Eq. (|14[) 
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(a) The critical load, p c , in a BA network. 
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FIG. 4. Numerical validation of the analytic results in BA (N = 512 nodes, m = 2) and ER (TV = 64 nodes, p = 0.1) networks. 
The continuous lines p c ,j and p c j represent Eqs. (|14[) and (|19[) . respectively. Numerical data, obtained by solving Eq. §5§ 
numerically, are represented by dots. (Color online.) 



and Eq. (|19|) on the same graph with the same absorp- 
tion level. 

The results are presented in Fig. @] It can be seen 
that the zeroth order approximation (|14p is valid only 
for small and large absorption levels. On the contrary, 
the first order approximation (|19|l fits the numerical data 
closely on the whole range of absorption levels. In order 
to emphasize the clear advantage of Eq. (|19p , the abso- 
lute difference between formulas (|14|) and (|19jl and the 
numerical data is also shown in Fig. 2) 



2. Bipartite graphs 

The nodes of a bipartite graph can be divided into two 
disjoint groups, Q\ and Q2, in such a way that the edges 
of the graph connect only nodes from different groups. 
Since trees are bipartite graphs, bipartite graphs are ex- 
tremely important from the practical point of view. 

Suppose that there are iVi nodes in Q\ and N2 nodes in 
Q2 , and denote the average degree and loading probability 



in Qi (i = 1,2) are dk = M/N{ and ■p i = J2keG,P k / N - 
respectively. 

The critical traffic load in bipartite graphs can be cal- 
culated similarly to the non-bipartite graphs. The cal- 
culation is based on the symmetry of the spectra of the 
transition matrix. Here, we only summarize the results 
for homogeneous loading, details of the calculation and 
the case of heterogeneous loading are presented in Ap- 
pendix [B] 

In the case of a bipartite graph we obtain that each 
partition defines a separate critical loading probability. 
For low or high absorption levels and homogeneous traffic 
load we obtain, analogously to Eq. (O, that 



(20) 



J_ - , 1 (1 - AQVji + (1 - M)/<*2 

(1) — 11 — "max: 



1 

72T 



1 (i - M )7rf 2 + (1 - M )M 

fi 2-n 



where dmax is the maximal degree in Q^. The critical 
loading probability of the whole network is the lesser of 
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the two: p c — min(pc 1 \pc 2 '*) ■ 

This result is very similar to the case of non-bipartite 
graphs. It can be seen that p c depends only on the 
absorption level, /x, and global properties of the graph, 
namely the mean di and the maximal degree dmax, which 
can be obtained from the degree sequences of Q\ and Q 2 
straightforwardly. 

In the case of intermediate absorption levels, the zeroth 
order approximation of critical traffic load, presented in 
([20| . has to be corrected. Similarly to non-bipartite 
graphs, the first order correction, obtained from the spec- 
tral decomposition of the transition matrix, is propor- 
tional to di/hi. The critical loading probabilities pi 1 ^ 
and p^ , including the first order corrections and corre- 
sponding to the two sub-components of a bipartite graph, 



1 

7 1 ! 



I 1 (1 - ^/dx + (1 
max< 



+ (1 - fi)di/hi 



di+ 




(21) 



-di+ 



Consequently, the critical loading probability of the 

whole network is p c = min (p^\p^j ■ 

In order to validate our results, we calculated the crit- 
ical load p c by inverting (JSJ) numerically. For validation, 
we used BA scale-free trees, grown by preferential at- 
tachment, and compared the numerical data with (|20[) 
and (|2l"j) with the same graph and absorption level. The 
results are shown in Fig. [SJ Eq. (l2"Tj) , which includes first 
order corrections, fits the numerical data closely not only 
at low and high absorption levels, like Eq. (f20|) . but also 
in the intermediate absorption range. 



D. Error estimation in the large graph limit 

We have seen that in the case of non-bipartite graphs 
the largest absolute eigenvalue is | Ai | = 1, whereas in case 
of bipartite graphs it is |Ai| = |Ajv| = 1. In the power 
series expansion (|15l) we considered the largest absolute 
eigenvalue precisely, and neglected the higher order terms 
for |A| < 1. In this section we discuss the validity of our 
approximation, and show numerically how the precision 
of our model depends on the graph properties. 

The error of the derived formulas depends on the mag- 
nitude of the higher order terms that we neglected in the 
power series expansion (I15[) . These terms can be bounded 
from above with the second largest absolute eigenvalue, 
that is max{|A2|, |Ai\r|} for non-bipartite graphs and A2 
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(a) The critical load, p c , as a fuction of the absorption level, fi. 
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(b) The absolute error of the critical load, Ap c , as a function 
of the absorption level, fi. 

FIG. 5. Numerical validation of the analytic results on a 
BA tree (N = 1024 nodes, m = 1). The continuous lines p c$ i 
andp Cj / represent Eqs. (|20[) and (|21[) . respectively. Numerical 
data, obtained by solving Eq. ([5| numerically, are represented 
by dots. (Color online.) 



for bipartite graphs. The smaller the second largest ab- 
solute eigenvalues are, the smaller the error of Eq. (fl"9j) 
and Eq. is. 

Although it is easy to manipulate the eigenvalues A& 
formally, it is difficult to see how the eigenvalues depend 
on the graph properties. In order to see more easily how 
the error depends on the graph properties, let us intro- 
duce the mixing time from the theory of random walk on 
graphs. 

The mixing time, r, is the expected time until a 
particle, performing random walk on a graph, reaches 
a stationary distribution. It can be shown that the 
mixing rate, the reciprocal of the mixing time, is pre- 
cisely equal to the second largest eigenvalue, i.e. 1/r = 
max{|A2|, I A at I } for non-bipartite graphs and 1/t = A2 
for bipartite graphs [29|. Therefore, it is plausible to 
conclude that for graphs which have small mixing rate, 
the error will be also small. 

There is no known formula in the literature on how the 



mixing rate depends on the graph parameters in general 
[30l | . However, our numerical experiments showed that 
the mixing rate depends strongly on the edge density, 
M/N = d/2. In particular the mixing rate decreases as 
the edge density increases. In Fig. [JJ we can see the rela- 
tive error of the derived formulas, Ap c /p c , as the function 
of various graph parameters. In the inset we can see the 
mixing rate as the function of the corresponding graph 
parameter. 

In the case of a BA network, for example, the edge 
density is M/N ~ to, where to is the number of edges 
connecting the new nodes to the graph in preferential 
attachment. This means that the mixing rate, 1/r, re- 
mains constant with small fluctuations if m is fixed, even 
if N — > 00. This phenomenon is the same in bipartite 
and non-bipartite graphs. For example, one obtains a 
BA scale- free tree, which is a bipartite graph, if to = 1, 
In this case the error of the derived formulas is also con- 
stant, and it will not decrease even in the thermodynamic 
limit. 
In Fig 



6 (a; 



and Fig. 

and to = 



6(b) 



we can see BA networks 
fixed, and N varied. These 



with m = 1 

cases correspond to a bipartite and a non-bipartite graph, 
respectively. We can see that in both cases both the 
relative error and the mixing rate tends to a fixed value 
as N — > 00. On the other hand, in Fig. 6(c) the size of 
the graph is fixed, and to is varied. We can see that the 
mixing rate increases as to increases and, at the same 
time, the relative tends to zero. 

The construction of ER networks is fundamentally dif- 
ferent from the BA graphs. In the case of ER networks, 
the edge density is M/N = p(N - l)/2 ~ pN/2. There- 
fore, the edge density increases, and the relative error 
tends to zero, if either p or N is increased, and the other 
parameter is fixed. Note that in case of p = 1 the ER 
network is a complete graph, in which case the derived 
formulas are exact. In contrast, the relative error of the 
derived formulas tends to a fixed value if pN is fixed. 

Numerical simulations carried out on ER networks are 
shown in Fig. [7J We can see that the simulation results 
confirm our assumption that the relative error of the de- 
rived formulas decreases if the density of the network 
increases, and remain fixed if the density is fixed. 



VI. CONCLUSION 

In this paper, we studied congestion phenomena in 
queuing networks. We analyzed how the critical point of 
the phase transition between free and congested phases is 
influenced by the topological properties of the network. 
In order to study the influence of the network structure 
on the traffic dynamics in an arbitrary network, we ne- 
glected congestion control mechanisms, and we modeled 
the particle transport by a simple Markovian random 
walk. 

In our model the critical traffic load in the network can 
be controlled by the absorption level of the particles. We 



0.5 p 
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0.15 -j 
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0.01 L 



A Pc/Pc 
A Pc,/Pc 
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N/10 3 
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(a) BA tree with m = 1 fixed and N varied. 



0.28 - Ap c( /p c 
A Pc,/Pc 
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N/10 3 
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N/10 3 



(b) BA network with m = 2 fixed and N varied. 



012 ■ A Pc,/Pc 




2 12 22 32 42 52 62 



52 62 



(c) BA network with N = 1024 fixed and m varied. 

FIG. 6. The relative error, Ap c /p c , of the derived formulas 
as a function of N and m in BA networks with absorption 
rate fi — 0.8. The black dots and red triangles show the 
error of the derived zeroth and first order approximations, 
respectively. The insets show the dependence of the mixing 
rate, 1/r on N and m. Data points were averaged over 32 
graph realization in each case. (Color online.) 



first derived Eq. (|14j) , the zeroth order approximation for 
the critical traffic load for low and high absorption levels. 
This result has been obtained by De Martino et al [l[ on 
a model with a congestion control mechanism included. 
Our result confirms their finding that at the critical point 
the details of the congestion control mechanism are less 
important in some cases. 
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(a) ER network with N = 1024 fixed and Per varied. 
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der formula, higher order corrections include not only the 
global properties of the degree sequence, i.e. the mean 
and maximum degree, but also the local information on 
the network structure. The improvement achieved by the 
higher order correction was validated by numerical sim- 
ulations. 

We also demonstrated that in the case of intermedi- 
ate absorption levels the structure of the network can 
have dramatic effects on the analytic behavior of the crit- 
ical point. We showed, in particular, that one must pay 
special attention when considering a bipartite graph, be- 
cause the spectra of bipartite graphs are symmetrical. 
We derived Eq. pTj) and showed that the critical point 
in a bipartite graph is the maximum of the critical points 
of its sub-components. 

Finally, we investigated the validity of our model. We 
showed that the higher order terms that were neglected 
during our calculations depend on the spectral gap, which 
can also be expressed by the mixing rate in the graph. 
We presented numerical arguments that the mixing ra- 
tio, that is the precision of our approximations, strongly 
depends on the edge density M/N. This empirical fact 
is well known in the mathematical community, but up to 
now, as far as we know, there is no rigorous proof of the 
phenomenon. We would like to examine this question in 
detail in our future work. 



(b) ER network with per = 0.5 fixed and N varied. 
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(c) ER with Nper = 49 fixed and N varied. 

FIG. 7. The relative error, Ap c /p c , of the derived formulas as 
a function of N and per in ER networks with absorption rate 
H = 0.8 (Figs. |7(a)1 and [7(b)] ) and fi = 0.5 (Fig. f7{cj\ . The 
black points and red triangles show the error of the derived 
zeroth and first order approximation, respectively. The insets 
show the dependence of the mixing rate, 1/r, on N and per- 
Data points were averaged over 32 graph realizations in each 
case. (Color online.) 
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Appendix A 

Let us suppose that E — P is not invertible, which 
is equivalent to the statement that P has a normal- 
ized eigenvector x with eigenvalue 1. In this case, us- 
ing the assumption that there is at least one node where 
X\ Pjj < 1 (e.i. there is absorption at least at one node), 



the following inequalities will hold: 



N 



N N 



N 



(Al) 



i— 1 j — 1 i=l 

which is a contradiction. 



In our paper we also showed that in the case of in- 
termediate absorption levels the zeroth order formula is 
not valid, and higher order corrections are needed. We 
derived Eq. (|19p which incorporates the first order cor- 
rections to Eq. (fbl| and improves the precision of the 
critical point considerably. In contrast to the zeroth or- 



Appendix B 

Using an indexing of the nodes that is suitable for the 
definition of bipartite graphs, every vector mentioned in 
the main text can be split into two parts, p = (p 1 ,p 2 ) T , 
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d = (d 1 ,d 2 ) T | = (£i,£ 2 ) T sucn that the N i i N 2) com- 
ponents of the first (second) part belong to nodes in Q\ 
(02)- The form of A and D are the following. 



, , B 

, B T 



D 



Di 
D 2 



(Bl) 



where B is an N\ x N 2 , Di is an N\ x N\ and D 2 is an 
N2 x N 2 matrix. The transition matrix has the following 
form: 



AD 





B T D 1 " 1 



BDn 



(B2) 



The structure of D 1/2 AD 1/2 is similar to A, the 
off-diagonal block matrices are D 1 1 ^ 2 BD 2 1 ' 2 and its 
transpose. A one line calculation shows, that if Ui — 
(ui,fe, M2,fc) T i s an eigenvector of D _1 ' 2 AD -1 ^ 2 with 
eigenvalue A&, then the vector (tii,fe, — w 2 .fc) T , is also 
an eigenvector with eigenvalue —Xk, so the spectra of 
D _1 ^ 2 AD _1 ' 2 is symmetric to the origin. One conse- 
quence of this symmetry is that if the number of nodes 



is even (odd), the kernel dimension of ' AD ' is 
also even (odd) . For the sake of simplicity, we study only 
bipartite graphs with an even number of nodes. For large 
networks, this has no serious consequence. Let us define 
the matrices 1^ and I^ 2 "*: 



r(i) 







u 2,k U 2 k 



(B3) 



This also can be split into three parts as in |HJ 
(i) , 2li 1) + (l- M )4 2) , 



A fc =0 



/x 2 — fx 
r(i) 



(2) 



v I^ + (l- M )A fc I^ 
Z - 1 l-(l-M) 2 A 2 fe ' 



(B9) 



A fe #0 



but here, summation runs over only the first half of the 
spectra. The power series expansion of the summands of 
the last term is 

2I« + 2(1 - fx)X k l[ 2) + 2(1 - / i) 2 A 2 I« + . . . (BIO) 

Dropping all the terms except the first gives 

2 (1 - /i) 2 L (1) + (1 - fx)l\ 2) , , 

E+-- w 1 V ?±J_ (Bll 

fi 2 — fi 

as the inverse of E - (1 - /i)D _1/2 AD~ 1/2 , so 

(E-D-^E+K 1 -^-^ 1 -^, (B12) 
fx 2 — fx 

where 3± and are the following matrices: 



J (i) = J_(d 1 lT 
1 M V d 2 \ T 2 



and 



(B13) 



and 



r(2) 











(B4) 



Then, the spectral decomposition of D 1//2 AD 1 ^ 2 is 

N N/2 



£A fe( i«+if) = $>4 : 



(2) 



(B5) 



fe=i 



fc=i 



It is easier to perform the calculation on the spectral 
decomposition of D _1 ' 2 AD -1 ' 2 instead of AD 1 to get 
an approximation of (E — P) _1 : 

(E - P)- 1 = D X / 2 (E - (1 - mJD-^AD- 1 / 2 )-^- 1 / 2 . 

(B6) 

The spectral decomposition of the factor in the middle 
of the r. h. s. is 



JY 



\ - u k u k 
fcl-Cl-AO* 



(B7) 



(2)_ 1 ( d X \ T 2 



Jl M \d 2 l{ 
This gives the values of $, l and £ 2 : 



t _„ , 1 (i-^) 2 PiM + (i-/j)p 2 M ^ 

Si — Pi H 7> a l> 

fi 2 — fx 

e „ 1 1 (i-^) 2 p 2 M + (i-/i)PiA ^ 

4 2 = P 2 + o 

/i 2 — fx 



(B14) 



(B15) 



Eq. [20] gives the final result for homogeneous loading. 

To get the first finite size correction, we have to use 
not only the first, but also the second term in lBlOl Using 
Eq. IB51 the correction term to the inverse of E — (1 — 
/i)D" 1/2 AD~ 1/2 appears to be (1 - /x)D~ 1/2 AD~ 1/2 - 

(2) 

2(1 — fj,)I\ , so the corrected formula for the inverse is 



E , 2 (l- Af ) 2 I< 1 > + (l- A ,)3l( 2 > 
fx 2 — fx 

+ (i- m )d- 1 / 2 ad- 1 / 2 , 



(B16) 



which, using the symmetry of the spectra, can be written 
in the form 



N/2 tW j- n ,,m t< 2 > 

h i-a-M)^ • 



(B8) 



and the inverse of E — P is 



(E " P) " E+ ^ —x (B17) 

+ (l- M )AD- 1 . 



^3t( 2 ) 
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The corrected values of £ x and £ 2 are 



€1 - Pi + - 

A* 



1 (1 - M) 2 P"i/di + (i - m) 3 p 2 M 



2 -/x 



€2 -P2 



+ (l-/i)BD~ 1 l 1 

1 (1 ~ ») 2 P 2 /d2 + (1 - M) 3 PiM 
/i 2 — /i 

(l-M)B T Dr 1 l 2 



d 2 + 



(B18) 



In the case of homogeneous loading probabilities, this 
leads to the appearance of the harmonic means: 

|m „ 1 1 (i - m) 2 M + (i - m) 3 M , 



p /i 2 — (i 

+ (1 - /J,)di/hi, 



6 



1 



1 (i~ M ) 2 /rf 2 + (i- M ) 3 /d 1 
/i 2 — fj, 

(1 - n)di/hi, 



(B19) 



and the individual critical loading probabilities in C/i and 
£2 are those that are in Eq. [2T] 



Appendix C 

The following algorithm calculates the order parame- 
ter. Suppose that p — pe, where e is a normalized vec- 
tor. If p is close to zero, the network is in the uncongested 
phase, the components of £(p) satisfy Eq.[5]and the num- 
ber of the uncongested nodes is equal to the number of 
the nodes in the system. If one increases p slowly, one 
finds that one or more nodes will surely have at least one 
waiting particle in their queues in the stationary regime, 



i. e. at least one becomes 1 at a certain value of p. 
Increasing p toward this value drives these nodes to the 
congested state, and the expected value of the growing 
length at the queues at these nodes in one time step is 
P e i + ^2ijPij — 1- On the other hand, these congested 
nodes send particles to their neighbors with rates equal 
to the corresponding element of P. These ideas suggest 
the following algorithm. 

1. Set P (0) = P, = e, = and s to a small 
positive number. 



2. In the fcth step, calculate the vector 
!E ( fe )( S ) = (E( fe )-pW)- 1 ( S y«+ Z 



Oh 



(CI) 



Starting from the last value at the (k — l)th step, 
increase s until one of the components of be- 
comes 1, or s becomes p. If the latter is the case, 
equate the components of £ to the corresponding 
components of x^ k \ If the former is true, set £j = 1 

(k) 

at the node where x\ is equal to one - this is the 
new congested node. Increase every component of 
2W with the corresponding element of P located 
in the column of the new congested node. Delete 
the rows and columns of the new congested node 

in E^, ~P^- k \ x^ and z 1 -^ . This gives matrices 
E (fc+i) ) p(fc+i) and vectors x (k+i) and z (k+i) for 

the (k + l)th step. 

3. If all components of £ are calculated, the order pa- 
rameter is 



v(p) 



(C2) 
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